MacWorld 1997 September/Macworld (1997-09).dmg/Serious Software/Cherwell Scientific Demos/pro Fit/pro Fit 5.0 demo (fpu).sea/pro Fit 5.0 demo (fpu)/Functions & Programs/•Gadgets/Linear interpolation of data

Text File | 1996-06-01 | 5KB | 168 lines

{
    This function accesses the current data set of a window and returns
 the linear interpolation between the data points. This function can be used
 to compare two data sets that were measured using slightly different steps
 for the x-coordinate. Based on the two data sets, it is possible to produce a
 table with measured data at selected x-coordinates, and so on.
 The function also offers the possibility of a multiplication factor and of x-
 and y-offsets.

 To use this function, choose "Add To Menu" from the Misc menu
 (or click the Add button).
 The function's name appears in the Func menu.

 Use it as you would use any other function. Click the "Choose Data" button
 in the parameter window to  set the data set used by the function.

 If you are looking at the function in the preview window, you must click
 "Redraw" whenever you change some data in its reference data set!
 The function has no way to know that you are changing its data under its nose!
}

function LinInterpol;

description 'returns the lin. interpolation of the current x- and y- data set in a given window.',
            '$BChoose Data$Offsets Δx,Δy are available. Data MUST be ordered';

defaults   a[1]:=0,active,'Δx     (x-offset)';
           a[2]:=0,active,'Δy     (y-offset)';
           a[3]:=1,active,'ampl   (mulitplies all y-vals.)';

var dataWindow,xc,yc,globMin,globMax,yMin,yMax:extended;
            iMax:integer; {last valid row number for data access}
            j1,j2,Lastj:integer;

{the following two functions are used to access the data window.
 They start looking in the data cell in row i. If in row i, both
 xCol and yCol cells contain a valid number they return this number,
 otherwise they increase i by one and look at the next data cell.
 This is repeated until a valid cell is found or all cells have been
 looked at. Both functions use the global variables xc and yc for
 the data columns, and iMax for the last row used.
 iMax is used (to spare time) instead of nrRows because it can
 be that many empty rows follow the data set
}

procedure Initialize;
begin
  dataWindow := 0;
  xc := 0; yc := 0;
end;

function xx(i:integer):extended;
begin
    if i<1 then xx:=-1/0 else
    begin
        while ((not dataok(i,xc)) or (not dataok(i,yc))) and (i<=iMax) do i:=i+1;
        if i<=iMax then xx:=data[i,xc] else xx:=1/0;
    end;
end;

function yy(i:integer):extended;
begin
    if i<1 then yy:=-1/0 else
    begin
        while ((not dataok(i,xc)) or (not dataok(i,yc))) and (i<=iMax) do i:=i+1;
        if i<=iMax then yy:=data[i,yc] else yy:=1/0;
    end;
end;

function check;
begin
    if pNumber=30000 then    {if the button was clicked}
    begin
        if dataWindow>0 then
        begin
            If GetWindowType(dataWindow)<>dataType then dataWindow:=0;
        end;
        Input('$WData Window',dataWindow,'$C$XX column',xc, '$C$YY column',yc);
        if (dataWindow>0) then SetCurrentWindow(dataWindow);
        SetDefaultCols(xc,yc,-1,-1);
    end;
    check:=ok;
end;

procedure first;
var counter:integer;
begin
if dataWindow>0 then
begin
    If GetWindowType(dataWindow)=dataType then
            SetCurrentWindow(dataWindow)
    else
    begin
            dataWindow:=0;
            SetCurrentWindow( FrontmostWindow(dataType) )
    end;
end;

if xc<=0 then xc:=xColumn;
if yc<=0 then yc:=yColumn;

 counter:=NrRows+1;
 repeat counter:=counter-1 until dataok(counter,xc) and dataok(counter,xc);
 GlobMax:=data[counter,xc];ymax:=data[counter,yc];iMax:=counter;j2:=iMax;
 GlobMin:=xx(1);ymin:=yy(1);j1:=1;Lastj:=j1;
 end;

procedure find(lastIndex:integer;target:extended);
{this procedure sets the global index-variables j1,j2 to point
 to the x-values bracketing the target value. The procedure first looks
 in the neighbourhood of the parameter lastindex. This means that
 if first looks at j1=lastindex, j2=lastindex+1. If this is not
 good, the procedure looks in the right direction (assuming that
 the data are sorted), making always wider steps until target is
 bracketed. Afterwards the copy j1,j2 is found with the bisection method.
}
 var    j:integer;
                 p:integer;
 begin
    p:=1;
    if xx(lastIndex)<=target then
    begin
        j2:=lastIndex+1;
        while xx(j2)<target do
        begin
            p:=p+p;
            j2:=j2+p;
            if j2>iMax then j2:=iMax;
        end;
        j1:=j2-p;
    end
    else
    begin
        j1:=lastIndex-1;
        while xx(j1)>target do
        begin
            p:=p+p;
            j1:=j1-p;
            if p<1 then p:=1;
        end;
        j2:=j1+p;
    end;

    while (j2-j1)>1 do
    begin
             j:=round((j2+j1)/2-0.4999999);
             if xx(j)>target then j2:=j else j1:=j;
     end;
     while (not dataok(j1,xc)) or (not dataok(j1,yc)) do j1:=j1+1;
     while (not dataok(j2,xc)) or (not dataok(j2,yc)) do j2:=j2+1;
end;

begin
 x:=x+a[1];

 if x<GlobMin then
         Find(Lastj,GlobMin)
 else if x> Globmax then
         Find(Lastj,Globmax)
 else Find(Lastj,x);

 Lastj:=j1;

 y:=data[j1,yc]+(x-data[j1,xc])/(data[j2,xc]-data[j1,xc])*(data[j2,yc]-data[j1,yc]);

 y:=y*a[3];
 y:=y+a[2];
end;